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ABSTRACT 

We present a ballistic description of the formation and propagation of the working 
surface of a relativistic jet. Using simple laws of conservation of mass and linear mo- 
mentum at the working surface, we obtain a full description of the working surface 
flow parametrised by the initial velocity and mass injection rate. This allows us to 
compute analytically the energy release at any time in the working surface. We com- 
pare this model with the results obtained numerically through a new hydrodynamical 
code applied to the propagation of a relativistic fluid in one dimension in order to 
test the limits of our study. Finally, we compare our analytical results with observed 
light curves of five long gamma ray bursts and show that our model is in very good 
agreement with observations using simple periodic variations of the injected velocity 
profiles. This simple method allows us to recover initial mass discharge and energy 
output ejected during the burst. 

Key words: hydrodynamics - relativity - galaxies:jets - quasars: general - gamma 
rays: bursts. 



1 INTRODUCTION 

Apparent superluminal knots observed along relativistic jets 
of quasars and micro-quasars are generally in terpreted as 
shock waves moving through the jet. It was first iReesI (jl966l ) 
who mathematically predicted the apparent superluminal 
motions observed in extragalactic jets, due to geometrical 
effects and relativistic velocities in the proper motion of 
knots inside them. The same author proposed that these 
observed knots were produce d by a vary ing velocity of the 
flow that moves along the jet (|Reeslll978h . Additionally, ob- 
servations of blazars have been carried through many years 
showing that the variability of the intensity and polarisation 
are most likely gener ated by a transverse shock propagating 
along the jet (see, e.g. Hagen-Thorn et a l. 2007; S pada et al 
l200ll: ISahavanathan fc Misral 1200*51 )" Recently, Ijamil et al 



1 20081 ) have proposed an internal shock model for micro- 
quasar jets in order to investigate particle acceleration and 
radiation production in these astrophysical objects. Another 
situation where internal shock waves appear is the fireball 
model for long Gamma Ray Bursts (GRBs). In such model, 
an effective mechanism for the generation of the observed 
gamma rays is the existence of internal shock waves along 
the associated jet which are caused by velocity v ariations 
of the outflow (|Rees fc Meszaroslll994lPiranll2005h . Despite 
the alternative explanations to the origin of internal shocks 
(e.g. iNaravan fc Kumarl 120081 ) the fireball model has been 



widely accepted and we show in this article how observa- 
tions match quite well with such picture. 

Detailed numerical models explaining the origin and the 
characteristic rad iation associated to i ntern a l shocks have 
been developed (iBlandford fc Konigll Il979l : iHugh es et al.l 
ll985l : lMarscherlll996l ). with special attention to the polar- 
isation of the observed radiation as a probe to the inter- 
nal shocks accelerating the material. In this article we ex- 
plore the formation of internal working surfaces propagating 
along a relativistic jet originated from the periodic varia- 
tion of the source velocity and/or mass injection. Looking 
at the one-dimensional propagation of mass particles, we 
present an extension of the non-relativistic model formu- 
lated by ICanto et~ai1 (120001 ). In that work, the formation 
and evolution of internal working sur faces is modelled for the 
ejecta correspond ing to stellar jets (|Raga fc Kofmanlll992l : 
iRaga et al.ll2004l ). This model has the advantage of provid- 
ing analytic expressions for the kinetic power radiated by 
the mass particles colliding inside the working surface. The 
extension of such analysis to the extreme relativistic regime 
is presented here. Assuming that the pressure gradients be- 
tween the fluid particles are negligible, and that radiation 
timescales are much shorter than the time it takes to form 
a particular working surface, we are able to recover ana- 
lytic expressions for the speed of the particles at both shock 
fronts and for the luminosity of the shocked gas. Our ana- 
lytic model is compared to our new numerical aztekas code 
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(www.aztekas.org) which completely solves the equations 
of hydrodynamics in a single dimension for the relativistic 
regime. In the light of our simulation, the limits for the va- 
lidity of our model are discussed. A final and most important 
test for our model is the comparison with gamma ray obser- 
vations from long GRBs. We present five cases in which the 
luminosity can be reproduced with our analysis and report 
on the inferred physical conditions of the ejecta that give 
way to the observed blast waves. 

The article is organised as follows. In section [2] the dy- 
namics of the setup are discussed and the analytical descrip- 
tion of the problem is presented. In section [3] we provide an 
example for the case of a constant mass discharge with an 
oscillating initial velocity. In section|4]we present the numer- 
ical method used to solve the equations in one dimension 
with planar symmetry. The comparison with the analytical 
model is presented through the example of section [3] As a 
test for the limits of our analytic model we also run sim- 
ulations where the pressure of the fluid is non-vanishing. 
In section [1] we compare the results of the numerical and 
analytic methods with each other and discuss the applica- 
tions of our model. Finally, in section [5] we use our simple 
analytical model to fit the light curves of five long GRBs. 



2 DYNAMICS OF RELATIVISTIC WORKING 
SURFACES 

The formation of shocks along the structure of a relativistic 
jet has been explained by several phenomena such as the 
presence of inhomogeneities in the surrounding media, the 
deviations and changes in the geometry of jets, and the tem- 
poral fluctuations in the parameters of the ejection (see e.g. 



Rees fc Meszaroslll994l;)Mendozall2000l ; lMendoza fc Longairl 
200 ll . 120021 ; I Jamil et al.ll2008l . and references therein). Here 



we are concerned with the last situation. When the speed 
of the emitted mass particles varies with time, a faster but 
later fluid parcel eventually hits an earlier but slower ejec- 
tion producing an initial discontinuity which gives rise to 
a working surface, i.e. a contact discontinuity surface, two 
shock waves, and two regions of shocked flow that must 
be at rest with respect to the contact, as the shocks re- 
cede from the contact surface in its frame of reference. 
The working surface travels along the jet with an aver- 
age speed u ws , as measured in the frame of the jet source. 
This picture for the formation of ra diation shock surfaces is 
known as the internal shock model (|Rees fc Meszaroslll994 
IDaigne fc Mochkovitchlll998l ). Although several extensions 
and particular aspects of the model have been presented in 
the l i terature (see e.g Panaitescu et al] 1 19991 ; ISpada et al.l 
l200ll ; ISahavanathan fc Misrall2005l )7 there are no simple an- 
alytical descriptions of this phenomenon. Here we present 
an analytical approximation to the formation of an internal 
working surface along a relativistic jet. We assume that the 
radiation timescales are small compa red to the charact eris- 
tic dynamical times of the problem (|Spada et al.ll200ll) . In 
consequence the pressure of the fluid is negligible and the 
collision is described ballistically. This assumption is valid if 
the flow within the jet is nearl y adiabatic and non-turbulent 
(|Sahavanathan fc Misrall2005l ). 

To follow the evolution of the working surfaces, we con- 
sider a source ejecting material in a preferred direction x 




Figure 1. When a fast velocity flow 2 moves over a slow velocity 
flow 1 , a working surface (represented with a curved line) moving 
with velocity v VB is generated as a result of the interaction. 



with a velocity v(t) and a mass ejection rate m(r), both 
dependent on time r as measured from the jet's source. 

Once the material has been ejected from th e source, we 
assum e it will flow in a free-stream way (see e.g. iRaga et al.l 
1990). This approximation is valid since the Mach number 
of the flow is large and, as mentioned before, we emphasise 
that the radiation processes which cool down the fluid occur 
in times scales shorter than other dynamical times associ- 
ated the to problem studied in this article (cf. ISpada et al.l 
l200ll ). The formation of a working surface is studied as the 
intersection of two different parcels of material ejected at 
times n = and T2 with flow velocities vi = v(n) and 
V2 = v(t2) — v-i + ctT2 respectively (see Figure [TJ, with 



:= (du/di 



If a > 0, the second parcel will eventually 



reach the first one. At time T2, the distance between both 
parcels is given by v 1 (j2 — T\) and thus the time t m , mea- 
sured in the reference frame of the central engine (i.e. the 
observer's frame), when both of them merge is given by 



t m =— W17 (vi) {l-v 
a 



a 



CtVlT2} , 

{l — 7 2 (i>i) av! r 2 } , 



where 7 (u) := 1/^/(1 — it 2 ) represents the Lorentz factor 
of the flow with velocity u, and we have assumed that the 
speed of light c = 1. The working surface is then formed at 
a distance 



d f = (t m + Ai)vi, 

from the central engine, where At = T2 — t\ . 

Following the non -relativistic formalism first proposed 
bv lCanto etahl (|2000l ). we assume that the working surface 
is thin and that there are no m ass losses within it (e.g. b y 
sideways ejection of material (see lFalle fc Ragalll993l . ll995f )). 
Since the flow is approximated as a free-streaming one, its 
velocity v(x,t) as a function of the x coordinate and time t 
is simply 



v{x,i) = v (t) = 



(1) 



This relation implies that the position a; ws of the work- 
ing surface from the downstream flow is given by 
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Vi(t - Tl), 



(2) 



and the one corresponding to the upstream flow takes the 
form 



X„ s = V 2 (t - T 2 ). 



(3) 



Consistent with the assumption that the flow is free- 
streaming, the velocity of the working surface is given by 
the velocity t> WB of its centr e of mass, which is determined 
by (jLandau fc Lifshit dll994h 



1 f T2 

VvB ^~M~J r y(v(s))m( s ) v ( s )ds, (4) 
where the weighted mass M 1 ejected between times n and 

T 2 is 



/ J (v(s)) rh(s) ds 

J n 



(5) 



Using formula Q for the velocity ti ws , it follows that 
the position of the working surface is given by 

1 f T ' 2 

x'ws = (t - t 2 )v vb + — / 7 {v(s)) m(s) v(s) (r 2 - s) ds. 

M l J ri 

(6) 

From equations ([2]) and ([3]) it follows that the position 
of the working surface as a function of the times n and t 2 
is also given by 



and the energy _B WS of the material inside the working sur- 
face, which is given by 



m7» 



(10) 



where the Lorentz factor 7 WS of the working surface material 
is such that 7^ s 2 = 1 — v 2 Js . 

If we assume now that the energy loss along the jet, 
E r = Eq — i5 ws is completely radiated away, then the lumi- 
nosity L := dE T /dt of the working surface is given by 



L = 



+ p-7ws72 {VwsV(T2) - 



■ 12 



iti(t 2 ) 
di/dr 2 

rh(n) dn f m 3 , . . 2 \ 

- d^ d^ \ 7ws + m; 7ws71 ( " wsW(ri) " " ws) _ 71 

(ii) 

where the Lorentz factors 7^2 := 1 — i' 2 (t"i,2) and, as we did 
before, we keep t 2 as the free parameter in the expansion. 
In consequence, the luminosity L in equation is found 
by writing down the expressions for n, v 2 , « ws and the 
derivatives dri/dr2 as well as di/dr2 as functions of the free 
parameter t 2 . 



3 A CONSTANT DISCHARGE FLOW 

As an example of our analytic description, let us consider 
the particular case of a constant discharge m and calculate 
the luminosity L that is obtained through simple oscillations 
of the particle emission speed. We assume that the injected 
velocity v has a periodic form given by 



V1V2 
v 2 — Vi 



(T2 



Tl 



(7) 



In the same manner, from the same set of equations, we 
calculate the time t as a function of ri and r 2 , giving 



T 2 V 2 — TiVl 
V 2 - Vl 



(8) 



For a given value of the position s ws , expressions ([2]), 
([3]) and ((6} establish a relation between times n and t 2 . 
Taking t 2 as a parameter, we can construct the position 
and velocity of the working surface as a function of t 2 and 
calculate the values of relevant quantities to the problem, 
such as the energy available on the moving working surface. 
Such a relation is one to one as long as the ejection speed 
v(t) increases monotonically. 

In order to calculate the amount of kinetic energy ra- 
diated away as the working surface moves, we take into ac- 
count the energy Eq the material had when it was ejected, 
which is well approximated bjQ 



£0 



f r 2 

I rh(s) 7 (v(s)) ds, 



(9) 



1 The term J -y(s)m(s) ds is subdominant in our problem as long 
as the variation of the velocity does not dramatically drop the 
speed to very small values with respect to the speed of light. 



v(t) = wo + 77 sinr. 



(12) 



in which the constant 77 <g; 1 . This type of oscillatory emis- 
sion speeds have been widely used for the description of 
internal shocks in both the Newtonian and the r elativistic 
cases (cf. lCanto et al.ll200d : IPanaitescu et al.ll 19991 ). For this 
example we choose a system of units in which rh — 1. In ad- 
dition we set the time unit so that the oscillation frequency 
is uj — 1. As a consequence of this assumption, the lumi- 
nosity L is such that its dimensions are the same as the 
dimensions of m, i.e. 



m , 



(13) 



and is thus dimensionless. 

If one assumes that the injected flow is highly relativis- 
tic, it is then possible to solve analytically equations 
(|11[) at 0(7 _1 ). However, the analytic expressions are long 
and cumbersome because, as opposed to the non-relativistic 
case, the Lorentz factor appears in this description ubiqui- 
tously. We have performed numerical integration of equa- 
tions Q through (|lip with a velocity profile given by equa- 
tion (|12l) . using the values vo = 0.9 and rj 2 = 0.09. The 
results are shown in Figure [2] and ar e very similar in shape 
to the ones obtained bv lCanto et al ] <|2000l ). 

The abrupt bump in the luminosity shows that the ki- 
netic energy must be radiated very effectively giving rise 
to emissions of high energy photons at various wavelengths 
depending on the strength of the shock. 
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0.94 




Figure 2. The figure shows the dimensionless velocity of the 
working surface t> ws and its dimensionless luminosity L as a func- 
tion of time for a period of flow oscillation with an input velocity 
given by equation \V2) and a constant discharge flow, under the 
assumption of a ballistic flow and conditions mentioned in sec- 
tion [3] All these quantities are measured in the rest frame of the 
jet source. 



The radiation mechanism and characterisation of the 
spectrum of particular objects will be the subject of 
a subsequent article. For a detailed study of the radi- 
ation mechanisms of internal shocks in GRBs see e.g. 



iDaigne fc Mochkovitchl (jl99ct ) and ISahavanathan fc Misral 
|2005T ) for the case of radiative knots along jets in AGNs. 



4 1-D NUMERICAL SOLUTION 

In order to compare our analytical calculations and the re- 
sults shown in Figure [2j we test the validity of our model 
and its approximations through a ID Relativistic Hydrody- 
namic (RHD) code developed by us for general purposes 
in astrophysical situations. The code uses a finite differ- 
ence method to solve the conti nuity, energy and momen- 
tum equations given by (s ee e.g. iBlandford fc McKeelll976l : 
iHidalgo fc Mendozalliooa ) 



d'yn 1 dynv 



dt 



dr 



dt 



e + v p 



i_ d_ 

r k Q r 



k 

r v 



0. 

p + e 



1 



= 0, 



d_ 

dt 



p + e 
'l-v 2 



+ 



k 2 P + e 
r v 



1 



dp 
dr 



= 0. 



(14) 
(15) 
(16) 



respectively, with k = 0, 1, 2 for planar, cylindrical and 
spherical geometry. Here n represents the particle number 
density, e = nm+e is the energy per unit volume of the fluid, 
m := 1 is the average mass per particle and e the internal 
energy per unit volume. Using a Bondi- Wheeler equation of 
state of the form 



= (K-l)e, 



(17) 



where re = 4/3 for an ultrarelativistic gas, we complete a 
system of equations for the functions n, p and v. The sys- 
tem can then be solved with suitable boundary and initial 
conditions for these quantities. 

The numerical method is based on a finite differ- 
ence scheme, with a fixed mesh and a fixed time step. 
At any time, th e solutions to equation s (114ll - (|16| l are 
obtained using a iMacCormack fc Paullavl (1 19721 ) integra- 
tion method. We have tested our code with traditional 
well known tests, such as the propagation of relativis- 
tic blast waves and shock tube problems described by 
iMartf fc MulTerl (|2003h . The code uses the GNU Scientific 
Library (www . gnu . org/ sof tware/gsl ) for many of its math- 
ematical computations, since this library is very well tested. 
More details on the features of this code will be published 
elsewhere. This code, named the aztekas code, is available 
on the Internet ( |www . aztekas . org| ) under a GNU General 
Public License (GPL), as described by the terms of the Free 
Software Foundation (www.gnu.org). 

In order to make a numerical comparison with the re- 
sults obtained in section [3] we use planar symmetry in equa- 
tions (|14p ~ (|16|) . with a null pressure model. At position x — 
of our domain, we inject at any time t matter with constant 
mass discharge m. This injected mass is assumed to have a 
velocity given by equation (|12p . For simplicity, we assummc 
that we have a single species of particles and so, we set a 
constant injected particle number per unit time flow h = 1, 
which in turn implies that the particle number density n at 
the point x = is given by 



n(t, x — 0) 



yu(t) 



yi-v 2 (t) 

v(t) ■ 



(18) 



In order to analyse a single flash of luminosity, such 
as the one described in Figure [2j the value of the veloc- 
ity is assumed constant after the time t — 2n, i.e. v(t > 
2-7T, x — 0) — 0.9. The initial conditions for the flow are 
chosen such that v(t = 0,x) = 0.9, p(t = 0,x) — 0.001 
and n(t = 0,a;) = y/l - v 2 (t = 0,x)/v(t = 0,x). The 
shock waves obtained by the varying velocity of the flow are 
captured numeric ally by introducing an artificial viscosity 
jBook et alJll975l ). Once the position of both shock waves 
are known, then at each time, the energy _E WS of the flow be- 
tween both shock waves (the working surface) is calculated 
as the sum nyAx, where the summation is done for each 
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numerical cell of width Ax that lies between both shocks. 
The energy of the input flow is calculated as in equation ((9} 
and the derivative that appears in the calculation of the 
luminosity L := d (Eo — E ws ) /At is performed numerically 
and softened usin g the flux-corrector method presented by 
iBook et al.l (|l975l ). Figure [3] presents the numerical result as 
compared with the analytic prescription described before. 
This shows that the analytical approximation is a good way 
to describe the dynamics and energetics of working surfaces 
formed by a varying input velocity. 

We have also performed calculations for two more cases 
in which the pressure p in the hydrodynamical equations has 
been written as (p, with £ = 0.01, 0.02. The pressure and 
the density are assumed to follow a polytropic relation of 
the form p oc n i ^ 3 , which is i n agreement w ith relation (|17p . 
and k = 4/3, as described bv lTooperl l|l965h . The initial and 
boundary conditions were not modified. However, the input 
energy Eo, given by equation ((9} is now 



E 



1 + (3 + v 2 ) 



1 dp 

n At 



(3 + v 2 ) 



Tdt, (19) 



according to equation l|15jl. The energy _E WS within the work- 
ing surface is calculated as the sum £ ws = nyAx + 
X^P7 (3 + v 2 ) Ax, where the summation is taken along all 
cells of width Ax of the domain, which lie between both 
shock waves. As it can be seen from the results presented in 
Figure [3j the peak of the luminosity is formed at the same 
time. However, the intensity of the pulse increases with an 
increasing £. The case f = 1, which corresponds to a full 
ultrarelativistic flow has not been drawn in Figure [3] since 
its luminosity peak has a much greater value. 
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Figure 3. The figure shows the dimensionless luminosity L of the 
working surface as a function of the dimensionless time t for a sin- 
gle period (t = 2ir) of oscillation with an input velocity given by 
equation H 12 ^ and a constant discharge flow. The solid line shows 
the analytical approximation under the assumption of zero pres- 
sure. From bottom to top, the dotted curves show the numerical 
computations made for the cases £ = 0,0.01,0.02. The reason 
as to why the luminosity doesn't go to zero at sufficiently large 
times (> 35) for the numerical solution is because the input ve- 
locity reaches a constant value = 0.9 after a period of oscillation, 
whereas the analytical approximation decays to a null value. 



5 ASTROPHYSICAL APPLICATIONS 

The model developed in the previous sections can be applied 
to shocks within jets emerging from AGNs, /i-QSRs and long 
GRBs. In order to see how this simple prescription can ac- 
count for some of the light curves of these objects, we select 
five long gamma-ray bursts and show that our luminosity 
function fits quite well their observed light curves. 

Our sample consist of five GRBs: GRB051111, 
GRB060206, GRB060904B,GRB070318 and 0GRB080413B 
with known redshi ft, observed by the BAT instrument on 
board the SWIFT jGehrelsl \20(m satellite (taken from the 
public database at ftp://legacy.gsfc.nasa.gov/swift), 
in the energy range from 15 - 150 KeV. In order to ob- 
tain the Flux F, the spectra taken at different sections of 
the light curve for each individual event were adjusted as a 
sample power law, with a normalisation N, and a photon 
ind ex a (for a more det ailed explanation of spectral analysis 
see lFirmani et al]l2008l . and references therein). 

In order to fit our model to the observations, we have 
made use of our null pressure analytical model with the same 
assumptions as the ones used in section[3] but with vo = 0.99 
and rj 2 — 0.009, and so the Lorentz factor of the injected flow 
varies from ~ 50 to ~ 500 in an oscillating sinusoidal way. 
To obtain the Flux F from the analytical approximation, we 
divide the analytical Luminosity L by 4nD 2 , where D is the 
luminosity distance, with cosmological parameters given by 
H = 71kms~ 1 Mpc -1 , fi ma ttor = 0.27 and S] vacmm = 0.73. 



Note that in order to get dimensional quantities, the dimen- 
sionless luminosity L has to be multiplied by mc 2 , with the 
unknown quantity rh. In other words, we will use the fact 
that that L D b s oc L, and so F a b B oc F. In the same manner, 
the dimensionless time t and the observer time t h B must be 
proportional to each other, i.e. t ^ s oc t. The proportionality 
factors are obtained by a linear regression analysis applied 
to both F's and t's separately. The results of these fits are 
shown in Figure [4] Note that the observed luminosity L ^ B 
represents an upper limit for the luminosity, since we have 
assummed that the efficiency factor e of converting injected 
kinetic energy to radiation has been tak en as one. In real- 
ity e < 1 and according to the results of lStern fc Poutanenl 
|200ct ) is close to one for Lorentz factors greater than 40. 

Complicated GRB light curve profiles may have to be 
adjusted by a mixture of a sum of sinusoidal variations of 
not only the velocity as shown in equation (I13|l . but also as 
a sum of periodic sinusoidal variations of the mass discharge 
rh, something outside the scope of this article. 



6 CONCLUSION 

We have constructed a full relativistic solution applied to 
the problem of a jet with varying periodic injection velocities 
and/or mass discharge. This was done by assuming that the 
working surfaces formed along the jet are ballistic. Under 
these circumstances we were able to obtain a full analytic 
description of their behaviour and with this, the luminosity 
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along the jet was calculated. We have also made numerical 
comparisons of our analytic approximations using a ID RHD 
code, the current status of the aztekas code. 

We prooved that the analytic approximations are in 
very good agreement with a full numerical solution under 
the assumption of null pressure gradients along the flow. 
Using the numerical code, the initial hydrodynamical quan- 
tities such as the particle discharge h, the particle number 
density n, the velocity v and the pressure p can be chosen 
to be intricated periodic functions of time. 

We have shown how to fit astronomical observations to 
five long GRBs using a simple regression analysis technique, 
assuming a constant discharge flow, a simple sinusoidal vari- 
ation of the velocity of the injected flow, and a ballistic flow 
approximation to the problem. It remains to test the analy- 
sis and match observations with more complicated shapes of 
luminosity curves. These could be produced by more com- 
plex flows like the sum of periodic sinusoidal functions for 
the injected velocity and mass discharge. Such possibilities 
are left for future research. 
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Figure 4. The figure shows the burst light curves (represented by dots with error bars at a la level) of five long GRBs, from top 
to bottom, left to right, GRB051111, GRB060206, GRB060904B, GRB070318, and GRB080413B observed with the BAT instrument 
onboard of the SWIFT satellite. The continuous curves on each graph are the best adjustments using the analytical model built in this 
article using a single period of oscillation for the variation of the velocity having a sinusoidal form given by equation 11121 (sec section 
[5] for more details). The fits give values of the isotropic luminosity given by (12.73, 141.8, 2.042, 1.771, 367.6) X 10 52 ergs s _1 for each 
burst respectively, which divided by the velocity of light squared imply mass injection discharges m ~ (1CT 1 - io- 2 )m q s- 1 . 



